library(rstan)
library(tidyverse)

# Set to the local directory of the replication data
setwd("~/Downloads/Replication_Data/")
n_cores <- 4

options(mc.cores = n_cores)

# Placement data from each year of the ANES
L2016 <- readRDS("Data/L2016.rds")
L2012 <- readRDS("Data/L2012.rds")
L2008 <- readRDS("Data/L2008.rds")
L2004 <- readRDS("Data/L2004.rds")
L2000 <- readRDS("Data/L2000.rds")
L1996 <- readRDS("Data/L1996.rds")
L1992 <- readRDS("Data/L1992.rds")
L1988 <- readRDS("Data/L1988.rds")
L1984 <- readRDS("Data/L1984.rds")
L1980 <- readRDS("Data/L1980.rds")
L1976 <- readRDS("Data/L1976.rds")
L1972 <- readRDS("Data/L1972.rds")

model_data_2016 <- list(n_items = length(unique(L2016$j)),
                        n_categories = length(unique(L2016$y)),
                        n_respondents = length(unique(L2016$k)),
                        n_obs = nrow(L2016),
                        s = as.numeric(table(L2016$j)),
                        k = L2016$k, j = L2016$j, y = L2016$y)

model_data_2012 <- list(n_items = length(unique(L2012$j)),
                        n_categories = length(unique(L2012$y)),
                        n_respondents = length(unique(L2012$k)),
                        n_obs = nrow(L2012),
                        s = as.numeric(table(L2012$j)),
                        k = L2012$k, j = L2012$j, y = L2012$y)

model_data_2008 <- list(n_items = length(unique(L2008$j)),
                        n_categories = length(unique(L2008$y)),
                        n_respondents = length(unique(L2008$k)),
                        n_obs = nrow(L2008),
                        s = as.numeric(table(L2008$j)),
                        k = L2008$k, j = L2008$j, y = L2008$y)

model_data_2004 <- list(n_items = length(unique(L2004$j)),
                        n_categories = length(unique(L2004$y)),
                        n_respondents = length(unique(L2004$k)),
                        n_obs = nrow(L2004),
                        s = as.numeric(table(L2004$j)),
                        k = L2004$k, j = L2004$j, y = L2004$y)

model_data_2000 <- list(n_items = length(unique(L2000$j)),
                        n_categories = length(unique(L2000$y)),
                        n_respondents = length(unique(L2000$k)),
                        n_obs = nrow(L2000),
                        s = as.numeric(table(L2000$j)),
                        k = L2000$k, j = L2000$j, y = L2000$y)

model_data_1996 <- list(n_items = length(unique(L1996$j)),
                        n_categories = length(unique(L1996$y)),
                        n_respondents = length(unique(L1996$k)),
                        n_obs = nrow(L1996),
                        s = as.numeric(table(L1996$j)),
                        k = L1996$k, j = L1996$j, y = L1996$y)

model_data_1992 <- list(n_items = length(unique(L1992$j)),
                        n_categories = length(unique(L1992$y)),
                        n_respondents = length(unique(L1992$k)),
                        n_obs = nrow(L1992),
                        s = as.numeric(table(L1992$j)),
                        k = L1992$k, j = L1992$j, y = L1992$y)

model_data_1988 <- list(n_items = length(unique(L1988$j)),
                        n_categories = length(unique(L1988$y)),
                        n_respondents = length(unique(L1988$k)),
                        n_obs = nrow(L1988),
                        s = as.numeric(table(L1988$j)),
                        k = L1988$k, j = L1988$j, y = L1988$y)

model_data_1984 <- list(n_items = length(unique(L1984$j)),
                        n_categories = length(unique(L1984$y)),
                        n_respondents = length(unique(L1984$k)),
                        n_obs = nrow(L1984),
                        s = as.numeric(table(L1984$j)),
                        k = L1984$k, j = L1984$j, y = L1984$y)

model_data_1980 <- list(n_items = length(unique(L1980$j)),
                        n_categories = length(unique(L1980$y)),
                        n_respondents = length(unique(L1980$k)),
                        n_obs = nrow(L1980),
                        s = as.numeric(table(L1980$j)),
                        k = L1980$k, j = L1980$j, y = L1980$y)

model_data_1976 <- list(n_items = length(unique(L1976$j)),
                        n_categories = length(unique(L1976$y)),
                        n_respondents = length(unique(L1976$k)),
                        n_obs = nrow(L1976),
                        s = as.numeric(table(L1976$j)),
                        k = L1976$k, j = L1976$j, y = L1976$y)

model_data_1972 <- list(n_items = length(unique(L1972$j)),
                        n_categories = length(unique(L1972$y)),
                        n_respondents = length(unique(L1972$k)),
                        n_obs = nrow(L1972),
                        s = as.numeric(table(L1972$j)),
                        k = L1972$k, j = L1972$j, y = L1972$y)

##### Compile to C

model_am <- stan_model(file = "AM_Sigma.stan")


##### Sample from posterior

# Parameters to exclude from posterior object to save memory
exclude_pars <- c("beta_raw", "alpha_raw", "zeta_raw")

# 2016
posterior_anes_2016 <- sampling(model_am, data = model_data_2016,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_2016, "Fitted_Models/posterior_anes_2016.rds")
rm(posterior_anes_2016); gc()

# 2012
posterior_anes_2012 <- sampling(model_am, data = model_data_2012,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_2012, "Fitted_Models/posterior_anes_2012.rds")
rm(posterior_anes_2012); gc()

# 2008
posterior_anes_2008 <- sampling(model_am, data = model_data_2008,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_2008, "Fitted_Models/posterior_anes_2008.rds")
rm(posterior_anes_2008); gc()

# 2004
posterior_anes_2004 <- sampling(model_am, data = model_data_2004,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_2004, "Fitted_Models/posterior_anes_2004.rds")
rm(posterior_anes_2004); gc()

# 2000
posterior_anes_2000 <- sampling(model_am, data = model_data_2000,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_2000, "Fitted_Models/posterior_anes_2000.rds")
rm(posterior_anes_2000); gc()

# 1996
posterior_anes_1996 <- sampling(model_am, data = model_data_1996,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1996, "Fitted_Models/posterior_anes_1996.rds")
rm(posterior_anes_1996); gc()

# 1992
posterior_anes_1992 <- sampling(model_am, data = model_data_1992,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1992, "Fitted_Models/posterior_anes_1992.rds")
rm(posterior_anes_1992); gc()

# 1988
posterior_anes_1988 <- sampling(model_am, data = model_data_1988,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1988, "Fitted_Models/posterior_anes_1988.rds")
rm(posterior_anes_1988); gc()

# 1984
posterior_anes_1984 <- sampling(model_am, data = model_data_1984,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1984, "Fitted_Models/posterior_anes_1984.rds")
rm(posterior_anes_1984); gc()

# 1980
posterior_anes_1980 <- sampling(model_am, data = model_data_1980,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1980, "Fitted_Models/posterior_anes_1980.rds")
rm(posterior_anes_1980); gc()

# 1976
posterior_anes_1976 <- sampling(model_am, data = model_data_1976,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1976, "Fitted_Models/posterior_anes_1976.rds")
rm(posterior_anes_1976); gc()

# 1972
posterior_anes_1972 <- sampling(model_am, data = model_data_1972,
                                warmup = 1000, iter = 5000, thin = 4,
                                control = list(adapt_delta = 0.999),
                                seed = 1, refresh = 1000,
                                include = FALSE, pars = exclude_pars,
                                open_progress = FALSE, chains = n_cores)
saveRDS(posterior_anes_1972, "Fitted_Models/posterior_anes_1972.rds")
rm(posterior_anes_1972); gc()

